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Summary. We review the theory of elliptic functions leading to SojiOTapes's for- 
mula for the sign function over the range e < |a;| < 1. We show how Gaufi' 
arithmetico-geometric mean allows us to evaluate elliptic functions cheaply, and 
thus to compute SojiOTapcB coefficients "on the fly" as a function of e. This in turn 
allows us to calculate the matrix functions sgn H, ^/H, and l/VH both quickly and 
accurately for any Hermitian matrix H whose spectrum lies in the specified range. 



1 Introduction 

The purpose of this paper is to provide a detailed account of how to compute the 
coefficients of SojiOTapes's optimal rational approximation to the sgn function. 
This is of considerable interest for lattice QCD because evaluation of the Neuberger 
overlap operator [1, 2, 3, 4] requires computation of the sgn function applied to a 
Hermitian matrix H. Numerical techniques for applying a rational approximation 
to a matrix are discussed in a companion paper [5], and in [6, 7]. 

In general, the computation of optimal (^eSbiniCB) rational approximations for a 
continuous function over a compact interval requires an iterative numerical algorithm 
[8, 9], but for the function sgnH (and the related functions and [5]) the 

coefficients of the optimal approximation are known in closed form in terms of Jacobi 
elliptic functions [10]. 

We give a fairly detailed summary of the theory of elliptic functions (§2) [11, 12] 
leading to the principal modular transformation of degree n (§2.7), which directly 
leads to SojiOTapes's formula (§3). Our approach closely follows that presented 
in [11]. 

We also explain how to evaluate the elliptic functions necessary to compute 
the 3ojiOTapcB coefficients (§4.4), explaining the use of the appropriate modular 
transformations (§2.7) and of Gaufi' arithmetico-geometric mean (§4.2), as well as 
providing explicit samples of code for the latter (§4.3). 
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2 Elliptic Functions 

2.1 Introduction 

There are two commonly encountered types of elliptic functions: Weierstrafi (§2.4) 
and Jacobi (§2.6) functions. In principle these are completely equivalent: indeed each 
may be expressed in terms of the other (§2.6); but in practice Weierstrafi functions 
are more elegant and natural for a theoretical discussion, and Jacobi functions are 
the more convenient for numerical use in most applications. 

Elliptic functions are doubly periodic complex analytic functions (§2.2); the com- 
bination of their periodicity and analyticity leads to very strong constraints on their 
structure, and these constraints arc most easily extracted by use of Liouville's theo- 
rem (§2.3). The constraints imply that for a fixed pair of periods uj and u' an elliptic 
function is uniquely determined, up to an overall constant factor, by the locations 
of its poles and zeros. In particular, this means that any elliptic function may be ex- 
panded in rational function or partial fraction form in terms of the Weierstrafi (§2.5) 
function with the same periods and its derivative. Furthermore, this in turn shows 
that all elliptic functions satisfy an addition theorem (§2.5) which allows us to write 
these expansions in terms of Weierstrafi functions with unshifted argument z (§2.5). 

There are many different choices of periods that lead to the same period lattice, 
and this representation theorem allows us to express them in terms of each other: 
such transformations are called modular transformations of degree one. We may also 
specify periods whose period lattice properly contains the original period lattice; and 
elliptic functions with these periods may be represented rationally in terms of the 
original ones. These (non-invertible) transformations are modular transformations 
of higher degree, and the set of all modular transformations form a semigroup that 
is generated by a few basic transformations (the .Jacobi real and imaginary trans- 
formations, and the principal transformation of degree n, for instance). The form of 
these modular transformations may be found using the representation theorem by 
matching the location of the poles and zeros of the functions, and fixing the overall 
constant at some suitable position. 

One of the periods of the Weierstrafi functions may be eliminated by rescaling 
the argument, and if we accept this trivial transformation then all elliptic functions 
may be expressed in terms of the Jacobi functions (§2.6) with a single parameter k. 

In order to evaluate the Jacobi functions for arbitrary argument and (real) pa- 
rameter we may first use the Jacobi real transformation to write them in terms of 
Jacobi functions whose parameter lies in the unit interval; then we may use the 
addition theorem to write them in terms of functions of purely real and imaginary 
arguments, and finally use the Jacobi imaginary transformation to rewrite the lat- 
ter in terms of functions with real arguments. This can all be done numerically or 
analytically, and is explained in detail for the case of interest in (§4). 

We are left with the problem of evaluating Jacobi functions with real argument 
and parameter in the unit interval. This may be done very efficiently by use of 
Gaufi' method of the arithmetico-geometric mean (§4.2). This makes use of the 
particular case of the principal modular transformation of degree 2, known as the 

Gaufi transformation to show that the mapping (a, 6) t-^ (^^{a + h),\fab^ iterates 
to a fixed point; for a suitable choice of the initial values the value of the fixed point 
gives us the value of the complete elliptic integral K, and with just a little more 
effort it can be induced to give us the values of all the Jacobi functions too (§4). 
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This procedure is sufficiently fast and accurate that the time taken to evaluate 

the coefficients of the SojioxapcB approximation for any reasonable values of the 
range specified by e and the degree n is negligible compared to the cost of applying 
the approximate operator sgn H to a vector. 

2.2 Periodic functions 

A function / : C — > C is periodic with period a; € C if f{z) = f{z + u)). Clearly, 
if a;i,a;2,-- - are periods of / then any linear combination of them with integer 
coefficients, UiUJi, is also a period; thus the periods form a Z-module. 

It is obvious that if / is a constant then this Z-module is dense in C, but the 
converse holds too, for if there is a sequence of periods aJi,a;2, • • • that converges to 
zero, then f'{z) = lim„_>oo [f{z + a;„) — f{z)] /uin = 0. It follows that every non- 
constant function must have a set of primitive periods, that is ones that are not 
sums of integer multiples of periods of smaller magnitude. Jacobi showed that if / is 
not constant it can have at most two primitive periods, and that these two periods 
cannot be colinear. 

2.3 Liouville's Theorem 

Prom here on we shall consider only doubly periodic meromorphic functions, which 
for historical reasons are called elliptic functions, whose non-colinear primitive pe- 
riods we shall call cj and co' . Consider the integral of such a function / around the 
parallelogram dP defined by its primitive periods. 



Substituting z' = z — u) in the second integral and z" = z — uj' in the third we have 



upon observing that the integrands identically vanish due to the periodicity of /. On 
the other hand, since / is meromorphic we can evaluate it in terms of its residues, 
and hence we find that the sum of the residues at all the poles of / in P is zero. Since 
the sum of the residues at all the poles of an elliptic function are zero an elliptic 
function cannot have less than two poles, taking multiplicity into account. 

Several useful corollaries follow immediately from this theorem. Consider the 
logarithmic derivative g{z) = [ln/(2:)]' = f{z)'/f{z) where / is any elliptic function 
which is not identically zero. We see immediately that g is holomorphic everywhere 
except at the discrete set {Qj} where / has a pole or a zero. Near these singularities 
/ has the Laurent expansion f(z) = Cj{z — C,jY^ +0 [{z — Cj)*^^^^) with Cj G C and 
Tj € Z, so the residue of g at C,j is rj. Applying the previous result to the function 
g instead of / we find that §gp dz g{z) = 2ni Y2j = 0, or in other words that the 
number of poles of / must equal the number of zeros of /, counting multiplicity in 
both cases. 

It follows immediately that there are no non-constant holomorphic elliptic func- 
tions: for if there was an analytic elliptic function / with no poles then f{z) — a 
could have no zeros either. 
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If we consider the function h{z) = zg{z) then we find 
(b dz h{z) = dz h{z) + dz h{z) + dz h{z) + dz h{z) 

JdP Jo J(jj Juj+ui' J u>' 

= f dz [h{z) - h{z + u')] + I dz [h{z + u:) - h{z)] 
Jo Jo 

= I dz [zg{z) - {z + u!')g{z)\ + I dz [{z + uj)g{z) - zg{z)] 
Jo Jo 

= —w' / dz g{z) + CO / dz g{z) 
Jo Jo 

= -Lo' {ln[/(w) - a] - ln[/(0) - a]} + co {ln[/(w') - a] - ln[/(0) - a]} 
= 2m{n'uj' + nuj), 

where n, n' € N are the number of times f{z) winds around the origin as z is taken 
along the straight hne from to a; or w'. On the other hand, Cauchy's theorem tells 
us that 

/ dzh{z)^I ^^Ifp. = 2nif]{a,-I3,), 
JdP Jds J[^) ^ 

where Uk and f3k are the locations of the poles and zeros respectively of /(«), again 
counting multiplicity. Consequently we have that J2T=i('^i' ~ Pk) = nuj + n'oj' , that 
is, the sum of the locations of the poles minus the sum of the location of the zeros 
of any elliptic function is zero modulo its periods. 



2.4 Weierstr£i6 elliptic functions 

The most elegant formalism for elliptic functions is due to Weierstrafi. A simple 

way to construct a doubly periodic function out of some analytic function / is to 
construct the double sum rn'e^ fi^~''nto—m'Lo). In order for this sum to converge 
uniformly it suffices that \f{z)\ < k/z^, so a simple choice is Q{z) = —2 m'ei^^~ 
muj — m'uj')'^. Clearly this function is doubly periodic, Q{z+u)) = Q{z+ijj') = Q{z), 
and odd, Q{—z) — —Q{z). 

The derivative of an elliptic function is clearly also an elliptic function, but in 
general the integral of an elliptic function is not an elliptic function. Indeed, if we 
define the Weierstrafi p function"' such that p' = Q we know that p{z+u!) = p{z) + c 
for any period lj. In this case we also know that p must be an even function, 
p{—z) = p{z), because Q is an odd function, and thus we have pi^to) = p{—^u!) + c 
by periodicity and p{^uj) = p{—^u!) by symmetry, and hence c = 0. We have thus 
shown that 

is an elliptic function. Its only singularities are a double pole at the origin and its 
periodic images. 

^ The name of the function is p, but I do not know what the name of the function is 
called; q.v., "Through the Looking-Glass, and what Alice found there," Chapter 
VIII, p. 306, footnote 8 [13]. 
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If we expand p in a Laurent series about the origin we obtain 

^^^^ z^^l^ 2^ (ma; + m'a;')2(^+i)"^2 + 20^28^ ' 

|™l + |™'l^o 

where the coefficients axe functions only of the periods 

— = 1 , _£3_ ^ \ - 1 

60 ^ (muj + m'uj'Y' 140 ^ (rruv + m'uj')^' ^ ' 



|m| + |m'|7^0 \m.\ + \m.'l^O 



l2 



From this we find p'{z) — —2z + ^aQ-zz + ifg:',z''' + • • •, and therefore [p'(2)] 
4«-«{l-iff22"-i53«« + ---} and [p{z)f = {l + ^jg^z* + ^53^" + • • •}• 
Putting these together we find 

[p'{z)f - 4[p{z)f + g2p{z) = -93 +Az'' + Bz^ + ---. 

The left-hand side is an elliptic function with periods uj and uj' whose only poles 
arc at the origin and its periodic images, the right-hand side has the value — at 
the origin, and thus by Liouville's theorem it must be a constant. We thus have 
[^'(0)]^ = 4[p(a:)]^ — g2p{z) — gs as the differential equation satisfied by p. Indeed, 
this equation allows us to express all the derivatives of p in terms of p and p'; for 
example 



p" = 6p2 - ifl2, p'" = 12pp', 

= 6(20p3 - 3g2P - 2g3), P^'^ = 18(20p2 - g^)p'. 



We can formally solve the differential equation for p to obtain the elliptic integral 
which is the functional inverse of p (for fixed periods u> and w'). 



Jo Jo 



dcp'io _ r 



\/4p(C)^ -fi'2p(C) -ff3 Jp{z) sj 4w3 _ - g-i 

It is useful to factor the cubic polynomial which occurs in the differential equa- 
tion, p'^(z) = 4(p — ei)(p — e'2)(<p — ea), where the symmetric polynomials of the 

roots satisfy 61-1-62-1-63 = 0, 6162 -I- 6263 -I- 6361 = —\gi, 616263 = \gz-, and 
e? + ei -I- 6^ = (ei -1-62-1- 63)^ - 2(6x62 -I- 6263 + 6361) = \gi. 

Since p' is an odd function we have — p'(ia;) = p'(— io;) = p'(ia;) = 0, and 
likewise p'(iw') = and p' (i(a; -|- w')) = 0. The values of p at the half-periods 
must be distinct, for if p(|w) = p(ia;') then the elliptic function p(«) — p(iw') 
would have a double zero at z = \lo' and at z = ^uj, which would violate Liouville's 
theorem. Since the p' vanishes at the half periods the differential equation implies 
that p{\ijj) = ei, p{\u)') = 62, p (i(ai -h w')) = 63, and that ei, 62, and 63 are all 
distinct. 

The solution of the corresponding differential equation with a generic quartic 

polynomial, y'"^ — a(y—ri){tj—r2){y—r'j,){y—r4,) with 7^ Vj , is easily found in terms 
of the Weierstrafi function by a conformal transformation. First one root is mapped 
to infinity by the transformation y — r4, + l/x, giving x'^ = —a{x — Pi){x — p2){x — 
pz) I pip2p?. with pj = l/(rj — Ti). Then the linear transformation x = + B with 
A = —Ap\p2pil a and B = (pi -|-p2-|-p3)/3 maps this to ^'^ = 4(5 — ei)({ — e2)(^ — 63) 
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where ej = {pj — B)/A. The solution is thus y = r4 + l/{Ap-\- B), where p has the 
periods impUcitly specified by the roots ej . 

It is not obvious that there exist periods to and u)' such that (72 and (73 are given 
by (1), nevertheless this is so (see [11] for a proof). 

A simple example of this is given by the Jacobi elliptic function snz, which 
is defined by z = f^'^^ dt [{1 — t^){l — k'^t^)] and hence satisfies the diff'er- 
ential equation {sn zf — (l + (snz)^) (l — fc'^(sn2)'^) together with the bound- 
ary condition snO = 0. We may move one of the roots to infinity by substitut- 
ing snz = 1 + l/x{z) and multiplying through by x{z)'^, giving x'^ = —2(1 — 
k'^) [x + ^][x — k/{l — k)][x + k/{l + k)]. The linear change of variable x{z) — 
— [125(2) -|- 1 — 5k^] I [6(1 — fc^)] then puts this into Weierstrafi's canonical form 
^'2 = 4(5 _ ei)(5 - e2)(5 - es) = - 9ii - 93 with the roots 

k'^ + 1 k'^±6k + l 
ei = ^— , ea±i= ; (3) 

and correspondingly 32 = A(fe'' -|- 14fe^ -|- 1), and 33 = -^{k^ + l)(fc^ + 6k + l)(fe^ - 
6k + 1). Clearly the Weierstrafi function p{z) with periods corresponding to the 
roots Cj is a solution to this equation. A more general solution may be written as 
^{z) = p{f{z)) for some analytic function /; for this to be a solution it must satisfy 
the differential equation, which requires that f'{z)^ = 1, so 5(«) = p(±« + ^) with 
A a suitable constant chosen to satisfy the boundary conditions. It turns out that 
the boundary values required for sn are satisfied by the choice ^{z) = p{z — K{k)), 

where K{k) = dt [(1 — i^)(l — fc^t^)] ^ is the complete elliptic integral. We shall 
later derive the expression for sn in terms of the Weierstrafi functions p and p' with 
the same argument and periods by a simpler method. 



The Weierstrafi ^-Function 

It is useful to consider integrals of p, even though these are not elliptic functions. If 
we define C' = ~Pj whose solution is 

c(.)4-/^.{p(.)-^} 

where the path of integration avoids all the singularities of p (i.e., the periodic 
images of the origin) except for the origin itself. The only singularities of C, arc a 
simple pole with unit residue at the origin and its periodic images. Furthermore, 1^ is 
an odd function. However, Q is not periodic: C,{z+lj) — (^{z)+ri where ?7 is a constant. 
Setting z = — and using the fact that C, is odd we obtain C,(—^u>) — Cik^) +''1 ~ 
—Ci^u)), or C{^Lo) = ^rj. If we integrate ^ around a period parallelogram P containing 
the origin we find a useful identity relating uj, w', r? and 77': 2m = §gpdzC,{z) = 

dz^ laz) - Ciz + a;')] + C^' dz [C(^ + a;) - C(«)] = C^' dz n - dz = 
770;' — rj'u). 



The Weierstrafi cr-Punction 



That was so much fun that we will do it again. Let (In a)' = a' /a = so 
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a{z) = zexp 



where again the integration path avoids all the singularities of ^ except the origin. 
(T is a holomorphic function having only simple zeros lying at the origin and its 

periodic images, and it is odd. To find the values of a on the period lattice wo 
integrate a'{z + uj)/a{z + uj) = (^{z + a;) = (^{z) + 7] = a'{z)/a{z) + r; to obtain 
\na{z + uj) = \na{z) + rjz + c, or cr{z + w) = c'e^'^(7{z). As usual we can find the 
constant c' by evaluating this expression ai z = —^u), (t{^iv) = —c'e~'^^^a{^(jj), 
giving c' = — ei**" and a{z + oj) = — e'''^^i"V(2;). 

2.5 Expansion of Elliptic Functions 

Every rational function R(z) can be expressed in two canonical forms, either in a 
fully factored representation which makes aJl the poles and zeros explicit, 

. ^ {z - bi){z - b2) ■ ■ ■ jz - b„) 
(z - ai)(z - a2) - ■ ■ (z - am) 

or in a partial fraction expansion which makes the leading "divergent" part of its 
Laurent expansion about its poles manifest. 

In these expressions hi are the zeros of R, ai its poles, c and yl^*' are constants, and 
S is a polynomial. It is perhaps most natural to think of E, the entire part of R, as 
the leading terms of its Laurent expansion about infinity. 

An arbitrary elliptic function / with periods cj and uj' may be expanded in two 
analogous ways in terms of Weierstrafi elliptic functions with the same periods. 

Multiplicative Form 

To obtain the first representation recall that ^^^lidj —bj) = (mod w, w'), so we 

can choose a set of poles and zeros (not necessarily in the fundamental parallelogram) 
whose sum is zero. For instance, we could just add the appropriate integer multiples 
of u! and to ai. We now construct the function 

^^^^ ^ a{z - bi)a{z - hi) ■ ■ ■ a{z - b„) 

a{z — ai)a{z — 02) • • ■ (j(z — an) ' 

which has the same zeros and poles as /. Furthermore, it is also an elliptic func- 
tion, since g{z + uj) = exp [^'7X]"=i('^j ~ ''j)] = sM- It follows that the ratio 
f{u) /g(u) is an elliptic function with no poles, as for each pole in the numerator there 
is a corresponding pole in the denominator, and for each zero in the denominator 
there is a corresponding zero in the numerator. Therefore, by Liouville's theorem, 
the ratio must be a constant f(u)/g{u) = C, so we have 

, , ^ a{z - bi)a{z - bj) ■ ■ ■ a{z - b„) 
a{z — ai)a(z — 02) • • • a-[z — an) 
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Additive Form 

For the second "partial fraction" representation let ai , . . . , a„ be the poles of / lying 
in some fundamental parallelogram. In this case, unlike the previous one, we ignore 
multiplicity and count each pole just once in this list. Further, let the leading terms of 
the Laurent expansion about z = akhe J2T^ii~^Y i''~~^y-^k^~^H^~ the func- 

flfe) then has exactly the same leading terms in its Laurent expansion. 

Summing this expression over all the poles, we obtain g{z) = X)fc=ifl'fc(^) = 
~ X]fc=i ^fe ^^C'' ^ o.k)- The sum of the terms with r > 1, being sums 

of the elliptic function p(z — Uk) and its derivatives, is an elliptic function. The 
sum of terms with r = 1, (p{z) = — X)"=i ■^k'^^Ci^ ~ Ofc)i behaves under translation 
by a period as ifi{z + lu) = ip(z) — '"/X]fe=i ^i"^ = H^{^)j where wc have used the 
corollary of Liouvillc's theorem that the sum of the residues at all the poles of the 
elliptic function / in a fundamental parallelogram is zero. It follows that the sum 
of r = 1 terms is an elliptic function also, so the difference f{z) — g{z) is an elliptic 
function with no singularities, and thus by Liouville's theorem is a constant C. We 
have thus obtain the expansion of an arbitrary elliptic function f{z) = C + g{z) = 

C — Y2^=i ^^'^~'^''C''^~^''('* ~ '^k), where the C functions have the same periods 

as / does. 



Addition Theorems 

Consider the elliptic function f{u) = p' (u) / (p(u) — p{v)); according to Liouvillc's 
theorem the denominator must have exactly two simple zeros, at which p{u) = p{v), 
within any fundamental parallelogram, p is an even function, p{—v) = p{v), so these 
zeros occur at it = ±v. At m = the function / has a simple pole, and the leading 
terms of the Laurent series for / about these three poles is (u— u)^^ + (u+v)~^ — 2/u. 
The "partial fraction" expansion of / is thus /(«) = C + (^(u — v) + (^{u + v) — 2(^{u), 
and since both / and C are odd functions we observe that C = 0. 

Adding this result, p'{u)/ {p{u) — p{v)) = C(m — v) + ^{u + v) — 2^(u), to the 
corresponding equation with u and v interchanged, —p'{v)/ {p{u) — p{v)) — —(^{u — 
«)+C(«+«)-2C(«), gives {p'{u) - p'{v)) I - p{v)) = 2C(«+v)-2C(w)-2C(^;). 

Rearranging this gives the addition theorem for zeta functions (^(u + w) = + 

av) + h - p'{v)) I - p(v)). 

The corresponding addition theorem for p is easily obtained by differentiating 
this relation 



-p{u + v) = -p{u) + 



1 p"{u)[p{u) - pjv)] - p'{n)[p'{u) - p'{v)] 



[piu) - p(vW 

and adding to it the same formula with u and v interchanged to obtain 

'-"{u) - p"iv)]lp{u) - piv)] - Ip'iu) - p'{v)]' 



-2p{u + v) = -p{u) - p{v) + i- 



[p{u) - p{v)Y 



Recalling that a consequence of the differential equation satisfied by p is the iden- 
tity (2), 2p" = 12p^ - 52, we have p"{u) - p"{v) = 6[p{uf - p{vf], and thus 

Pin +v) = -p{u) - p{v) + 1 
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Differentiating this addition theorem for p gives the addition theorem for p'. 
Since higher derivatives of p can be expressed in terms of p and p' there is no need 
to repeat this construction again. 

Representation of Elliptic Functions in terms of p and p' 

Consider the "partial fraction" expansion of an arbitrary elliptic function / in terms 

of zcta functions and their derivatives, f{z) — C + X]fc=i A^i^~^\''^~^^ {z — at)- 

This expresses f{z) as a linear combination of zcta functions 1^(2; — at), which are 
not elliptic functions, and their derivatives (^^^\z — ak) = — p^^^''\z — at) which are. 
We may now use the addition theorems to write this in terms of the zeta function 
(^{z) and its derivatives C^^\z) = —p^^~^\z) of the unshifted argument z. 

For the r — 1 terms the zeta function addition theorem gives us X]fc=i ^k'^Ci^ " 
ak) = X)fc=i -^k'^Ci^) + Ri {p{^)j p'{^))j were we use the notation Ri{x, y) to denote 
a rational function of x and y; i.e., an element of the field C{x,y). The coeffi- 
cients in R\ depend transcendentally on a^, of course. This expression simplifies to 
just the rational function i?i(p, p') on recalling that, as we have previously shown, 

Using the addition theorems for p and p' all the terms for r > 1 may be ex- 
pressed in the form ELi ^'^'^C^'^'H^ " Ofc) = " ELi 4''"' V'''^'^^ " afe) = 
Rr {p{z). p' [z)). We have thus shown that / = i?(p, p'). In fact, since the differen- 
tial equation for p expresses p'^ as a polynomial in p we can simplify this to the 
form / = Re{p) + Ro{p)p' ■ A simple corollary is that if / is an even function then 
/ = Rc(p) and if it is odd then / = Ro{p)p' . 

A corollary of this result is that any two elliptic functions with the same periods 
are algebraic functions of each other. If / and g arc two such functions then / = 
Ri{p) + R2{p)p', g = Ri.{p) + Ri{p)p' , and p'^ = 4p^ - g2p - gr, these equations 
immediately give three polynomial relations between the values /, g, p, and p' , and 
we may eliminate the last two to obtain a polynomial equation F(f,g) = 0. To be 
concrete, suppose Ri{z) = Pi{z) /Qi{z) with Pi,Qi G C[2], then we have 

Pi(/, p, p') = Qi{p)Q2{p)f - Pi{p)Q2{p) - P2{p)Qi{p)p' = 0, 
P2{g,p,p') = Q3{p)Q4{p)g - P3{p)Qa{p) - Pa{p)Q3{p)p' = 0, 

■PsCp, p') = p'^ - 4p^ -I- gip -I- 33 = 0; 
we may then construct the resultants^ 

P4(/,p) = Res(Pi(/,p,p'),J^?(p,p')) =0, 
P5(3,p) = Rcs{P2{g,p,p'),Pz{p,p')) =0, 
F{f,g) = Res(P4(/,p),P5(ff,p)) =0. 

A corollary of this corollary is obtained by letting g = /', which tells us that every 
elliptic function satisfies a first order differential equation of the form F{f, /') = 
withF€C[/,/']. 

^ In practice Grobner basis methods might be preferred. 
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A second metacoroUary is obtained by considering g{u) = f{u + v), for which 

we deduce that there is a polynomial C{v)[f(u)][f{u + v)] 3 F = 0, where C{v) 
is the space of complex- valued transcendental functions of v. On the other hand, 
interchanging u and v we observe that F £ C{u)[f{v)][f{u + v)] too. The coefficients 
of F are therefore both polynomials in /(w) with coefficients which are functions 
of V, and polynomials in f{v) with coefficients which are functions of u. It therefore 
follows that the coefficients must be polynomials in f{u) and f{v) with constant 
coefficients, F £ C[f{u),f{v),f{u + v)]. In other words, every elliptic equation has 
an algebraic addition theorem. 



2.6 Jacobi Elliptic Functions 

We shall now consider the Jacobi elliptic function sn implicitly defined hy z = 
dt [{1 — t^){l — k^t^)] ^. This cannot be anything new — it must be express- 
ible rationally in terms of the Weierstrafi functions p and p' with the same periods. 

The integrand of the integral defining sn has a two-sheeted Riemann surface 
with four branch points. The values of z for which sn(2, k) = s for any particular 
value s € C are specified by the integral; we immediately see that there two such 
values, corresponding to which sheet of the integrand we end up on, plus arbitrary 
integer multiples of the two periods w and w'. These periods correspond to the non- 
contractible loops that encircle any pair of the branch points. There are only two 
independent homotopically non-trivial loops because the contour which encircles all 
four branch points is contractible through the point at infinity (this is a regular 
point of the integrand, as may be seen by changing variable to 1/z). 

We may choose the first period to correspond to a contour C which contains the 
branch points at z = ±1. We find 



dt 



c^il-t^){l-kH^) 



dt 



V(l-i2)(l-fe2t2) 



= AK{k), 



taking into account the fact that the integrand changes sign as we go round each 
branch point onto the other sheet of the square root. Likewise, we may choose the 
second period to correspond to a contour C" enclosing the branch points at z = 1 
and z = 1/k; this gives 



dt 



l/k 



1/fc 



dt 



V'(l-i2)(l-fe2t2) 



If we change variable to s = -^^(1 — k'^t'^)/{l — k'^) we find that 



i/fe 



dt 



^{l-t^){l-kH2) Jo V(l-s2)(l-fc'2s2) 



ds 



iK(k'), 



where we define k' = \fY^^W. We thus have shown that the second period w' = 
2iK{k') is also expressible as a complete elliptic integral. 

The locations of the poles of sn are also easily found. Consider the integral 
J-^i^ dt [(1 — t^){l — k^t'^)]^^ , by the change of variable s — 1/kt wo see that it is 

equal to ds (fes^)"^ [(l - l//c^s^) (l - ~* = K{k). We therefore have that 
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Jo ^il-t^){l-kH^) 

and thus sn has a pole at 2K(k) + iK{k'). 

We mentioned that there are always two locations within the fundamental par- 
allelogram at which sn(2;, k) = s. One of these locations corresponds to a contour Ci 
which goes from t = on the principal sheet (the positive value of the square root 
in the integrand) to i = s on the same sheet, while the other goes from t = on the 
principal sheet to t = s on the second sheet. This latter contour is homotopic to one 
which goes from t = on the principal sheet to t = on the second sheet and then 
follows Ci but on the second sheet. If the value of the first integral is z, then the value 
of the second is 2K{k) — z, thus establishing the identity sn(2:, k) = sn(2if (fc) — z, k). 

Since the integrand is an even function of t the integral is an odd function of sn, 
from which we immediately see that 8n{z^ k) = — sn(— a, k). 

We summarise these results by giving some of the values of sn within the fun- 
damental parallelogram defined by = AK and lo' — 2iK' : 



z 





K 


2K 


3A" 


iK' 


K + iK' 


2K + iK' 


SA" + iK' 


sn(«, k) 





1 





-1 


00 


1/k 


— 00 


-1/k 



where we have used the notation K = K(k) and K' = K{k'). 



rl /'l/fc roc 

/ + / +/ 

Jo Jl Jl/k 



dt 



^{l-t^){l-kH^) 



: 2K{k) + iK{k'), 



Representation of sn in terms of p and p' 

From this knowledge of the periods, zeros, and poles of sn we can express it in terms 
of Weierstrafi elliptic functions. Prom (3) we know that the periods uj = 4K, uj' = 
2iK' , and lo+lj' = 4:K+2iK' correspond to the roots ei, 62, and 63; that is p {^uj) — 
ei, p (jw') = 62, and p {^{cl' + 1^')) = 63. Since sn{z, k) is an odd function of z it must 
be expressible as R{p{z)) p' {z) where 7i is a rational function; since it has simple 
poles in the fundamental parallelogram only aX, z = iw' = iK' and z = 1(01-1-0;') = 
2K -h iK' it must be of the form sn{z^ k) = R {p{z)) p'{z)/ [{p{z) — 62) {p{z) — ea)] 
with R a rational function.'' The WcierstraB function has a double pole at the origin, 
its derivative has a triple pole, and the Jacobi elliptic function sn has a simple zero, 
so we can deduce that R (p{z)) must be regular and non-zero at the origin, and hence 
.R is just a constant. As sn(«, k) = z + 0{z^) near the origin this constant is easily 
determined by considering the residues of the poles in the Weierstrafi functions, and 
we obtain the interesting identity 8n.{z, k) = —\p'{z)l [(p(a) — 62) — 63)]. 

Representation of sn^ in terms of p 

We can use the same technique to express sn(«, fe)^ in terms of Weierstrafi elliptic 
functions. The differential equation satisfied by 5(2) = sn(2;, A;)^ is s'^ = 4s(l — 

s)(l — fc'^s), which is reduced to Weierstrafi canonical form (^'^ = 4("'* — (/2C ^ .93 with 
^2 = K/c* -I- fc^ -I- 1) and gz = ^(2A:^ - 1)(A:^ - 2)(fc^ -|- 1) by the linear substitution 
s = (C + 1(1 + fe^)) /k"^ ■ The roots of this cubic form are ei = — i(l -|- A;^), esii = 

Remember that the logarithmic derivative of a function din f{z)/dz = f'{z)/f{z) 
always has a simple pole at each pole and zero of /. 
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i(2 ± fe^). The general solution of this equation is <^{z) = p{izz + c) where c is a 
constant, and since sn{z,k)'^ has a double pole at z = iK{k') wc have sn{z,k)'^ — 
[p (z — iK{k')) + ^(1 + fc^)] /k^ . This can be simplified using the addition formula 
for Weierstrafi elliptic functions to give** sn(2:,fc)^ = 1/ [p{z) — ei]. Of course, we 
could have seen this immediately by noting that sn(z, fe)^ is an even elliptic function 
with periods 2K{k) and 2iK{k') corresponding to the roots e;, and therefore must 
bo a rational function of p{z). Since it has a double pole at z = iK{k') and a 
double zero at « = in whose neighbourhood sn(«, k)^ = z^ + 0{z'^) the preceding 
expression is uniquely determined. 

A useful corollary of this result is that we can express the Weierstrafi function 
p{z) with periods 2K(k) and 2iK'{k) rationally in terms of sn{z, k)^, namely p{z) = 
sn{z, k)^^ + ei, and thus any even elliptic function with these periods may be written 
as a rational function of sn(z,fc)^. 



Addition Formula for Jacobi Elliptic Functions 

We may derive the explicit addition formula for Jacobi elliptic functions using a 
method introduced by Euler. Consider the functions si = sn(it, fc), S2 = sn{v,k) 
where we shall hold u + v = c constant. The differential equations for si and S2 are 
s'l = (1 — s?)(l — fc^sf), s'2 = (1 — si)(l — fc^si), where we have used a prime to 
indicate differentiation with respect to u and noted that v' = —1. Multiplying the 
equations by and respectively and subtracting them gives W{si,S2) ■ (S1S2)' = 
(siS2 — S2s'i)(siS2 + S2S1) = (,s? — si)(l — fc^sfsi), where we have introduced the 

Wronskian W{si,S2) = det ( ) . If we differentiate the differential equations for 

\ Si S2 J 

si and S2 wc obtain s'/ — —{1 + k'^)si + 2k'^ Si, s'2 — — (1 + fc'^)s2 + 2fe^S2; subtracting 
these equations gives W' = (siS2 — S2s'i)' = S1S2 — S2s'l — — 2fe^siS2(si — si) = 

may combine the expressions wc have derived for W 
and W to obtain (InM^)' = W /W = {l-k^sisl)' /{l-k^s{sl) = (ln(l - A;^s?si))'. 
Upon integration this yields an explicit expression for the Wronskian, W = C(l — 
fe'^Sisi) where C is a constant, by which wc mean that it docs not depend upon 
u. The constant does depend on the value oi c = u + v, and it may be found by 
evaluating formula at w = 0. 

To do so it is convenient to introduce two other Jacobi elliptic functions 
cn(u, fc) = \/l — sn(u, where cn(0, fc) = 1; and dn(u, fc) = ^/l^^Wsn(^uJe)^ , 
where dn(0, k) — 1. In terms of those functions wc may write sn' u = cnudnw, fur- 
thermore they satisfy the identities (snw)^-l-(cn-u)^ = 1 and (fc^ snu)^-l-(dnM)^ = 1, 
and differentiating these identities yields cn'u = — snwdn u, dn' u = —k^ snwcnw. 

We may now write C = W{snu,snv)/ \l — {ksnusnv}^^ = [snMcntidn« + 
snwcnwdnu]/ [l — (ksnusnv)'^], remembering that v' = —1. Setting v = gives 
C = snu = sn c, and thus we have the desired addition formula sn(u + v) = 
(sn-u cnw dn« + sn u cn u dn m) / (l — (fc sn« sn u)^) . 

* The Weierstrafi functions in this expression implicitly correspond to the roots ei, 
whereas as those in the previous expression for sn(2,fc) corresponded to the a. 
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2.7 Transformations of Elliptic Functions 

So far we have studied the dependence of clhptic functions on their argument for 
fixed values of the periods. Although the Woierstrafi function appear to depend 
on two arbitrary complex periods uj and uj' they really only depend on the ratio 
r = w'/w. If we rewrite the identity p{z) = p{z + u!) = p{z + u)') in terms of the new 
variable = z/cu wo have p{u>C,) ~ P {^{C + 1)) = P ("^(C + '''))• Viewed as a function 
of C we have an elliptic function with periods 1 and r, and as we have shown this 
is expressible rationally in terms of the Weierstrafi function and its derivative with 
these periods. 

Another observation is that there are many choices of periods uj and a;' which 
generate the same period lattice. Indeed, if we choose periods uj = auj + fSui' , uj' = 



between elliptic functions with these periods called a first degree transformation. 
Jacobi Imaginary Transformation 

Jacobi's imaginary transformation, or the second principal' first degree transforma- 
tion, corresponds to the interchange of periods u>' = —Cj and u> = ui' . We start with 
the function sn{z,k)^ which has periods u) = 2K and u' = 2iK' , and consider the 
function sn(2;/M, A)^ with periods uj = 2ML and lo' = 2iML' (with L = K{X) and 
L' = K{\') as usual). For suitable M and A we have ML = iK' and iML' = —K, 
corresponding to the desired interchange of periods. 

Since sn(2;/M, A)^ is an even function whose period lattice is the same as that of 
sn{z, kY it must be expressible as a rational function of sn(z, fe)^, and this rational 
function may be found by matching the location of poles and zeros. sn{z/M, A)^ has 
a double zero at z/M — and a double pole at z/M = iL' . This latter condition 
may be written as z = iML' = —K, or z = K upon using the periodicity conditions 
to map the pole into the fundamental parallelogram. Thus 



The constant A may be found by evaluating both sides of this equation a,t z = 
iK': on the left sn{iK' /M, X)^ = sn(L, A)^ = 1, whereas on the right wo have A 
because sn{z, k) ^ oo as z ^ iK' . We thus have A = 1. 

The value of A is found by evaluating both sides a,t z = —K + iK': on the 
left sn{{-K + iK')/M,xf = sn{iL' + L,Xf = 1/A^, and on the ri ght we have 
A/{1 - P) since sii{-K + iK' , kf = l/fc^. Wo thus have A = Vl - fc^ = k' . 

From these values for A and A we may easily find M, as iK' = ML = 
MK(X) = MK{k') = MK' gives M = i. We may therefore write the Jacobi 
imaginary transformation as sn{—iz, k')^ = sn{z, k)^ / (sn(«, k)^ ~ l) > or equiva- 
lently sn(i2, k') — i sn(z, k) / cn(2;, fc), where we have made use of the fact that sn^ is 
an even function, and chosen the sign of the square root according to the definition 
of cn and the fact that sn(2;, k) = z -\- 0{z'^). 

The first principal first degree transformation may be derived similarly. See [11] 
for details. 



70; + (5a;' with det 



a 13 
7 5 



= 1 then this will be the case. This induces a relation 
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Principal Transformation of Degree n 

We can also choose periods iD and uj' whose period lattice has the original one as 
a sublattice, for instance we may choose uj — cj and w' = co'/n where n € N. 
Elliptic functions with these periods must be rationally expressible in terms of the 
WeierstraB elliptic functions with the original periods, although the inverse may not 
be true. This relationship is called a transformation of degree n. 

Let us construct such an elliptic function with periods 4:K and 2iK' /n, where 
K = K{k) and K' = K{k') with k'^ + k''^ = 1. We may do this by taking sn{z,k) 
and scaling z by some factor 1/M and choosing a new parameter A. We are thus 
led to consider the function sn(z/M, A), whose periods with respect to z/M are 
4L = 4:K{\) and 2iL' = 2iK{\'), with A^ + A'^ = 1. Viewed as a function of z it 
has periods 4LM and 2iL' M, and M and A are thus fixed by the conditions that 
LM = K and L' M = K' /n. The ratio f{z) = sn{z/M,\)/ sn{z,k) must therefore 
be an even function of z with periods 2K and 2K'\ 



f{±z + 2mK + 2im'K'^ - ^ ' 



sn(±z + 2rnK + 2irn'K', k) 
^ sn{±^+2mL + 2im'nL',X) ^ ±(-)"' sn A) ^ 
sn{±z + 2mK + 2im'K',k) ±(-)'" sn(z, fc) 

for m,m' G Z; and hence f{z) must be a rational function of sn(2:, k)^ . 

Within its fundamental parallelogram the numerator of /(z) has simple zeros 
at « = m2iL'M = vn2iK' /n for m = 0, 1, . . . , n — 1 and simple poles for m = 
i,|,...,n— i; whereas its denominator has a simple zero at a: = and a simple pole 
at z = iK' . Hence, if n is even f(z) has simple zeros for m = 1, 2, . . . , in — 1, in + 
1, . . . , n — 1, a double zero for rn — in, and simple poles for m — i,|,...,n — i; 
whereas if n is odd then m = 1, 2, . . . , n— 1 give simple zeros and m = i, |, . . . , |n — 
1, in + 1, . . . , n — i simple poles. Therefore there are 2[inJ zeros and poles, and it 
is easy to see that they always come in pairs such that the zeros occur for an{z, k) = 
±sn{2iK'm/n,k) and the poles for sn(2;, fc) = ±sn{2iK' {m — ^)/n,k) with m = 
1, . . . , L|nJ. We thus see that the rational representation is 

/ 2 ,\ 1 sn(z,fc)^ 

f(^\ = Im'^J - _ it sn(2ijf'm/n,fc)^ , 

■'^ sn{z,k) M 11 1 ^ ' 

Bn(2z/f'(m-^)/n,fc)^ 

where the overall factor is determined by considering the behaviour near z = 0. 

The value of the quantity M may be determined by evaluating this expression 
at the half period K where sn{K, k) = 1 and sa.{K/M, A) = sn(L, A) = 1, so 

M=nr 



sn(2iK' m/n,k)^ 



sn(2iif'(m-l)/n,fc) 



The value of the parameter A is found by evaluating the identity at z = K + iK' /n, 
where sn (i^+i^, a) = sn (f + ii^, a) = sn(L + zL', A) = i. 
It will prove useful to write the identity in parametric form 
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with ^ = sn{z,k), Cm = sn(2iK'm/n, k)~^ and c'^ = sn{2iK' {m — ^)/n,k)~^ . 
This emphasises the fact that sa.{z/M, A) is a rational function of sn(^;, k) of degree 
(2LinJ + l,2LinJ). 



3 SojiOTapes's Problem 

SoJiOTapeB's fourth problem is to find the best uniform rational approximation to 
sgna; = — i9(— .x) over the interval [—1, — e] U [e, 1]. This is easily done using the 
identity (5) derived in the preceding section. 

We note that the function ^ = sn(«, k) with fc < 1 is real and increases monoton- 
ically in [0, 1] for z € [0, if], where as before we define K = K{k) to be a complete 
elliptic integral. Similarly we observe that sn(2:, k) is real and increases monotoni- 
cally in [1, 1/k] iox z = K + iy with y £ [0, A"] and K' = K{k'), k^ + k''^ = 1. On the 
other hand, sn(2;/M, A) has the same real period 2K as sn{z,k) and has an imagi- 
nary period 2iK' /n which divides that of sn{z, k) exactly n times. This means that 
sn(2/M, A) also increases monotonically in [0, 1] for z € [0,if], and then oscillates 
in [1, 1/A] for z^K + iy with y G [0, K'\. 

In order to produce an approximation of the required type we just need to rescale 
both the argument ^ so it ranges between —1 and 1 rather than — 1/fc and 1/k, and 
the function so that it oscillates symmetrically about 1 for ^ £ [1, l/fe] rather than 
between 1 and 1/A. We thus obtain 

2 X k^ -c 

X m=l 

with k = e. On the domain [—1, — e] U [e, 1] the error e{x) = R{x) — sgn(a;) satisfies 
\e(x)\ < A = j-^, or in other words \\R — sgn ||oo = A. Furthermore, the error 
alternates 4|_inJ + 2 times between the extreme values of ±Zi, so by MeSbinieB's 
theorem on optimal rational approximation R is the best rational approximation of 
degree (2[inJ + 1, 2[inJ). In fact we observe that R is deficient, as its denominator 
is of degree one lower than this; this must be so as we are approximating an odd 
function. Indeed, we may note that R'{x) = [1— A^) / R{x) is also an optimal rational 
approximation. 



4 Numerical Evaluation of Elliptic Fiinctions 

We wish to consider GauB' arithmetico-gcometric mean as it provides a good means 
of evaluating Jacobi elliptic functions numerically. 

4.1 Gaufi Transformation 

GauB considered the transformation that divides the second period of an elliptic 
function by two, u}[ = uji and 0)2 = ^u)2. This is a special case of the principal 
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transformation of nth degree on the second period considered before (4) with n = 2, 
hence 

sn(2,fc)'^ 
sn(iK' .k)-^ 



sn{z, ) 
M 



1 



sn(2,fc)'^ 
sn(iiif',fc)2 



with the parameter A corresponding to periods L = K/M and L' = K' /2M. Using 
Jacobi's imaginary transformation (the second principal first degree transformation, 
with Lo'i = —u!2 and W2 = wi), 



sn(i2;, k) = i 



sn{z, k') 
cn{z, k') ' 



(7) 



we get 



sn{z, k) 
M 



1 + 



1 + 



sn(z,fc)^ cn(K' fc')^ 
sn(z,fc)2 cn(iX',fc')2 



[i(iif',fe')2 

Since sn(J!",fe') = 1, cn(Ji",fc') = 0,® sn(iJi",fc') = l/x/T+T, and cn(iK',fc') = 
+ k), wo obtain sn(2/Af, A) =sn{z,k) [l + k sn{z, kf]~^ /M . 
To determine M weset z = K: 1 = sn{K/M, A) = sn(A', A;) [l + fe sn{K, kf] ~)m 
= [1 + k]~'^/M, hence M = 1/(1 + fc). To determine A we set m = + iK' /2: 



, K iK' , 
sn — + — TT , A = 
' M 2M' ' 



sn(A' + 
M[l + A;sn(i<:+ iiJ£:',fc)2]" 



Now, from the addition formula'^ sn(w -\- v) = (sn u cn v dn v + snt;cnwdnw)/[l — 
(fcsnwsnv)^], we deduce that 



sn[K + 



iK' 



sn cn if^ dn if^ + sn if^ cn /s: dn 



1 - {ksnKsn'-^y 



cn if^ dn if^ 



Furthermore sn{^iK' ,k) = isn{^K' ,k')/ cn(^K' ,k') = i/Vk, and correspondingly 
cn{^iK',k) = Y^(l + k)/k, and dn{iiK',k) = VI + fe, giving sn(Js: + \iK\k) = 
l/\/fc. We thus find 1/A = sn(L + iL' , A) = \/2Ms/k or A = 2M\/fc = 2v%/(l + k). 
Combining these results we obtain an explicit expression for Gaufi' transformation 



sn (1 + k)z. 



2Vk 
1 + k 



(1 + k) sn{z,k) 
H-fesn(z,fe)2 ■ 



^ Let X = sn(ijf, fc), then by the addition formula for Jacobi elliptic functions 
sn{K,k) = 1 = 2a;v'l-a:Vl-fc^a;V(l - k'^x*)- Hence (1 - fc^a;'')^ = 4x^(1 - 
a;^)(l - fc^a;^), so fc*2,-* - 4fc2a;'' + 2(2 + fc ^)a:^ - 42;^ + 1 = or, with a = - 1, 
[z^-(l-fe^)]^ =0. Thus 2 = ±Vl - fc2 = ±fe',ora; = 1/x/T^T'. SinceO < x < 1 
we must choose the positive sign, so sn {^K{k), k) = 1/Vl + k'. 

^ We shall suppress the parameter k when it is the same for all the functions 
occurring in an expression. 
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4.2 Arithmetico-Geometric Mean 



Let an,b„ a M. with a„ > hn > 0, and define their arithmetic and geometric means 
to be ttn+i = 5(0™ + bn), hn+i = Vo-nb-n- Since these are means we easily see that 
o„ > a„+i > bn and a„ > b„+i > b„; furthermore a^+i — 6^+1 = + 2o„b„ + 

b'i) — anbn = 3(0^ — 2anbn + 6n) = ^(ffln — &n)^ > 0, SO «„ > On+1 > &„+l > bn- 

Thus the sequence converges to the arithrnetico-geometric mean a^o ~ boo- 

If we choose a„ and bn such that k = {an — bn)/{an + bn), e.g., an = 1 + fc and 
6„ = 1 — fe, then 



1 + 



(On + bn) + (On — bn) _ On 



fc2 ^ / an - fen \ ^ (an + &n) - 4a„b„ ^ ^ 

V an + 6n / (a„ + &n)^ a^+i ' 



+ (l + fc)2 Vl + fc 

^ _ I " (an +&n) - (an - bn) V 



{an + 6n) + (an — &n) 

If we define Sn = sn ^(1 + k)z, j and Sn+i = sn(a, fe) then Gauf5' transformation 
tells us that 

_ (1 + k)Sn+l _ anSn + 1 _ 2anSn+l 

Sn 



1 + ^4+1 an+l [l+ (f^^) S^+i] (an + 6n) + (an - 6n)4+l ' 



On the other hand 



^(l-t^)[ 
and these two integrals may be rewritten as 



f '_!f _ ttn + l f 

J I r 7 ^2 \ T an J 

l{l-t^) 



r dt 1 r dt 

J 7(1 - t2)(l - A;2t2) ~ 1 + kJ ■ 



IT+W 



(If 



{l-t^)[l 



Therefore the quantity 



an+i 



dt r 

Jo ^/a-t^)\al^Al-t^) + bl^,m Jo 



dt 



^{1 - t^)[al+-,{l - t^) + bl^-^t^] Jo ^J{l-t^)[al{l-f^) + bnf^] 
is invariant under the transformation {ambn,Sn) i— > (on+i, &n+i, Sn+i), and thus 
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This implies that Soo = sin ^ j = sin(aoo2) with our previous choice of a„ = 

l + fe,fe„ = l — fe=> a„+i = l,bn+i = VI — /c^- We may thus compute Sn+i = 
sn{z, k) = f{z, 1, Vl-A;2) for < fe < 1 where 

f sin(o«) if o = 6 

" I with ^ ^ / ^ , V^) if a ^ 

Furthermore, if we take z — K{k) then s„+i = 1 and s„ = 2ari/[(a„ + 6„) + (on — 
&„)] = 1; thus Soo = sin(aooJ^) = 1, so UooK = tt/2 or K{k) = 7r/2aoo- 



4.3 Computer Implementation 

An implementation of this method is shown in Figures 1 and 2. 

The function cirithgeom recursively evaluates the function / defined above. One 
subtlety is the stopping criterion, which has to be chosen carefully to guarantee 
that the recursion will terminate (which docs not happen if the simpler criterion 
b==a is used instead) and which ensures that the solution is as accurate as possible 
whatever floating point precision FLOAT is specified. Another subtlety is how the 
value of the arithmetico-gcometric mean *agm is returned from the innermost level 
of the recursion. Ideally, we would like this value to be bound to an automatic 
variable in the calling procedure sncndnK rather than passed as an argument, thus 
avoiding copying its address for every level of recursion (as is done in here) or 
copying its value for every level if it were explicitly returned as a value. Unfortunately 
this is impossible, since the C programming language does not allow us to have 
nested procedures. The reason we have written it in the present form is so that 
the code is thread-safe: if we made agm a static global variable then two threads 
simultaneously invoking sncndnK might interfere with each other's value. The virtue 
of this approach is only slightly tarnished by the fact that the global variable pb used 
in the convergence tost is likewise not thread-safe. The envelope routine sncndnK is 
almost trivial, except that care is needed to get the sign of cn(2;, k) correct. 



4.4 Evaluation of SojiOTapeB Coefficients 

The arithmetico-geometric mean lets us evaluate Jacobi elliptic functions for real 

arguments z and real parameters < fe < 1. For complex arguments we can use the 
addition formula to evaluate sn{x + iy, k) in terms of sn(a;, k) and sn(iy, k), and the 
latter case with an imaginary argument may be rewritten in terms of real arguments 
using Jacobi's imaginary transformation. We can either use these transformations 
to evaluate elliptic functions of complex argument numerically, or to transform alge- 
braically the quantities we wish to evaluate into explicitly real form. Here wc shall 
follow the latter approach, as it is more efficient to apply the transformations once 
analytically. 

SoJiOTapes's formula (6) is 



T 2 2 

1 + 1 feM ■'■■I- fe2 _ c' 

A m=l 
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#include <math.h> 
#define ONE ((FLOAT) 1) 
ttdefine TWO ((FLOAT) 2) 
#define HALF (ONE/TWO) 

static void sncndnK (FLOAT z, FLOAT k, FLOAT* sn, FLOAT* cn, 
FLOAT* dn, FLOAT* K) { 

FLOAT agm; 

int sgn ; 

*sn = arithgeom(z, ONE, sqrt(ONE - k*k) , feagm) ; 
*K = M_PI / (TWO * agm) ; 

sgn = ((int) (fabs(z) / *K)) 4; /* sgn = 0, 1, 2, 3 */ 
sgn "= sgn » 1; /* (sgn & 1) = , 1 , 1 , */ 
sgn = 1 - ((sgn & 1) « 1); /* sgn = 1, -1, -1, 1 */ 
*cn = ((FLOAT) sgn) * sqrt(ONE - *sn * *sn) ; 
*dn = sqrt(ONE - k*k* *sn * *sn) ; 

} 



Fig. 1. The procedure sncndnK computes sn(2;,fe), cn(«,fe), dn{z,k), and K{k). It 
is essentially a wrapper for the routine arithgeom shown in Figure 2. The sign of 

cn(z, k) is defined to be —1 if K{k) < z < 3K{k) and +1 otherwise, and this sign is 
computed by some quite unnecessarily obfuscated bit manipulations. 

static FLOAT arithgeom (FLOAT z, FLOAT a, FLOAT b, FLOAT* agm) i 
static FLOAT pb = -ONE; 
FLOAT xi; 

if (b <= pb) I pb = -ONE; *agm = a; return sin(z * a) ; } 
pb = b; 

xi = arithgeom(z, HALF*(a+b), sqrt(a*b), agm); 
return 2*a*xi / ((a+b) + (a-b)*xi*xi) ; 

} 



Fig. 2. Recursive implementation of GauB' arithmetico-geometric mean, which is the 
kernel of the method used to compute the Jacobi elliptic functions with parameter 
k where < fc < 1. The function returns a value related to sn{z,k'), and also sets 
the value of *agm to the arithmetico-geometric mean. This value is simply related 
to complete elliptic function K{k') and also determines the sign of cn(^, fe'). The 
algorithm is deemed to have converged when b ceases to increase: this works whatever 
floating point precision FLOAT is specified. 
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with Cm = sn{2iK'm/n,k) ^ and = 8n.{2iK'(m — \)/n,k) ^. We may evaluate 
the coefficients Cm and by using Jacobi's imaginary transformation (7), 



cn{2K'm/n,k') 



sn{2K'm/n,k') 



i(2X'(r 



i)/n,k' 



sn(2K'(;m- ^)/n,k') 



We also know that M = nit"! (1 " c™)/(l - d^), and 1/A = (^7^^) nLi"i (1 " 
CmP)/(l — c'mi^) with f = sn{K + iK'/n,k). We may use the addition formula 
to express the Jacobi elliptic functions of complex argument in terms of ones with 
purely real or imaginary arguments, so 



sn( K+ , 

n 



cn dn 



sn /s: cn dn ^ + sn ^ cn if dn X 



1- (fcsn/sTsnif^) 



iK' 



l_(fcsnif^) dn^ 



These may be converted to expressions involving only real arguments by the use of 
Jacobi's imaginary transformation (7), 



sn I , k 
n 



cn I , k 
n 



dn ( , k 
n 



\ 



1 + 



\ 



1 + 



dnfi£l,A.' 



giving the simple result ^ = 1/ dn(/sr'/n, fe'). 
Putting these results together we have 



with 



Cm 



2g'(m-l) ^, 



Li"J 



2 1 -TT Cm / 1 - \ ^ ^ 1 - A 
1 + 1/AA: ^\ c'^ U-c™; ' 1 + A' 

where Zi is the maximum error of the approximation. 
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